##############################
# STATE REACH AND DEVELOPMENT // PSRM // Replication
#
# Descriptive plots and maps of state reach in space and time
# Sourced from "replication.R"
#
##############################




# Figure 1: Plot descriptive maps

# ... plot years
plot.yrs <- c(seq(1966, 2006, by = 20), 2016)


# ... plot globals
zlim = log(1+c(0, 160))
leg.ticks <- c(0, 1, 3, 7, 20, 50, 150)
plot.size <- list(large = c(6,6), small = c(2,2))
colpal <- inferno(n = 256, direction = -1)

# ... cshapes 2016
cshp.2016 <- cshp(date = as.Date("2016-01-01"))
cshp.2016 <- cshp.2016[cshp.2016$gwcode %in% afro.cow,]

# ... plot
for(t in c("regcap", "natcap")){
  
  # ... load stack
  stack <- raster::stack(file.path("data", paste0("statereach_", t, "_", plot.yrs, ".tif")))
  
  # ... log transform
  stack <- log(1+stack)
  
  # ... plot
  for(y in c(1:length(plot.yrs))){
    ## plot size?
    if(plot.yrs[y] == 2016){
      s = 1
    } else {
      s = 2
    }
    
    ## Plot
    png(paste0(fig.path,"fig1_descr_", t, "_", plot.yrs[y], "_", names(plot.size)[s], ".png"),
        width = plot.size[[s]][1], height = plot.size[[s]][2], units = "in", res = 300)
    par(mar = c(0,0,0,0))
    print(raster::plot(stack[[y]], useRaster = T, interpolate = F, axes=FALSE, bty="n", box=FALSE,
                       col = colpal, zlim = zlim,
                       legend = FALSE #, main = plot.yrs[y]
    ))
    if(names(plot.size)[s] == "large"){
      # ... country borders
      plot(as_Spatial(cshp.2016), add = T, border = "white")
      
      # ... legend
      r.range <- c(minValue(stack[[y]]), maxValue(stack[[y]]))
      plot(stack[[y]], legend.only=TRUE,  col = colpal, zlim = zlim,
           legend.width=1, 
           axis.args=list(at=log(1+leg.ticks),
                          labels=leg.ticks),
           legend.args=list(text='Travel time (hrs)', side=4, font=2, line=2.5, cex=.8))
    }
    dev.off()
  }
  
  # Cleanup
  rm(stack)
}

# Figure 3: National Capital ###################
# also contained in Figure A9 Appendix. 

# ... Load country-level aggregates
dist.cow.yrs <- read.csv(file.path("data", "dist_natcap_bycountry.csv"))

# ... Aggregate
dist.agg.df <- aggregate.data.frame(list(natcap.dist = dist.cow.yrs$natcap.dist * dist.cow.yrs$pop,
                                         pop = dist.cow.yrs$pop),
                                    by = dist.cow.yrs[,c("year"), drop = F], 
                                    FUN = sum, na.rm = T)
dist.agg.df$natcap.dist <- dist.agg.df$natcap.dist / dist.agg.df$pop


# ... plot
png(paste0(fig.path, "fig3_natcapdist.png"), width = 3.5, height = 3, units = "in", res = 300)
g <- ggplot(dist.agg.df[dist.agg.df$year >= 1960,], 
            aes(x = year, y = natcap.dist)) +
  geom_line(data = dist.cow.yrs[dist.cow.yrs$year >= 1960,], 
            aes(group = cowcode), col = "grey", alpha = .5)  + 
  geom_line() + 
  xlab("Year") + ylab("Average distance to\nnational capital (hrs; road-travel)") +
  theme_minimal()
print(g)
dev.off()


# ... numbers cited in text
writeLines(as.character(round_any(dist.agg.df$natcap.dist[dist.agg.df$year == 1965] - 
                                    dist.agg.df$natcap.dist[dist.agg.df$year == 2015],.1)) ,
           paste0(num.path,"incr_natcap_abs.tex"))
writeLines(as.character(round_any(dist.agg.df$natcap.dist[dist.agg.df$year == 1965],.1)) ,
           paste0(num.path,"incr_natcap_1965.tex"))
writeLines(as.character(round_any(dist.agg.df$natcap.dist[dist.agg.df$year == 2015],.1)) ,
           paste0(num.path,"incr_natcap_2015.tex"))
writeLines(as.character(round_any(100 * (dist.agg.df$natcap.dist[dist.agg.df$year == 1965] - 
                                           dist.agg.df$natcap.dist[dist.agg.df$year == 2015]) /
                                    dist.agg.df$natcap.dist[dist.agg.df$year == 1965],.1)) ,
           paste0(num.path,"incr_natcap_perc.tex"))


# Figure 3 and A9: Regional Capital ###################
# also contained in Appendix Figure A7  

# ... Load country-level aggregates
regcap.cow <- read.csv(file.path("data", "dist_regcap_bycountry.csv"))

# ... Aggregate
dist.agg.df <- aggregate.data.frame(list(regcap.dist = regcap.cow$regcap.dist * regcap.cow$pop,
                                         pop = regcap.cow$pop),
                                    by = regcap.cow[,c("year"), drop = F], 
                                    FUN = sum, na.rm = T)
dist.agg.df$regcap.dist <- dist.agg.df$regcap.dist / dist.agg.df$pop

# ... plot
png(paste0(fig.path, "fig3_regcapdist.png"), width = 3.5, height = 3, units = "in", res = 300)
g <- ggplot(dist.agg.df[dist.agg.df$year >= 1960,], 
            aes(x = year, y = regcap.dist)) +
  geom_line(data = regcap.cow[regcap.cow$year >= 1960,], 
            aes(group = cowcode), col = "grey", alpha = .5)  + 
  geom_line() + 
  xlab("Year") + ylab("Average distance to\nregional capital (hrs; road-travel)") +
  theme_minimal()
print(g)
dev.off()


# ... numbers cited in text
writeLines(as.character(round_any(dist.agg.df$regcap.dist[dist.agg.df$year == 1965] - 
                                    dist.agg.df$regcap.dist[dist.agg.df$year == 2015],.1)) ,
           paste0(num.path,"incr_regcap_abs.tex"))
writeLines(as.character(round_any(dist.agg.df$regcap.dist[dist.agg.df$year == 1965],.1)) ,
           paste0(num.path,"incr_regcap_1965.tex"))
writeLines(as.character(round_any(dist.agg.df$regcap.dist[dist.agg.df$year == 2015],.1)) ,
           paste0(num.path,"incr_regcap_2015.tex"))
writeLines(as.character(round_any(100 * (dist.agg.df$regcap.dist[dist.agg.df$year == 1965] - 
                                           dist.agg.df$regcap.dist[dist.agg.df$year == 2015]) /
                                    dist.agg.df$regcap.dist[dist.agg.df$year == 1965],.1)) ,
           paste0(num.path,"incr_regcap_perc.tex"))



# FIGURE A8: Plot within-country variation in increase

# ... prepare data

# ... --- load country borders
country <- cshp(date = as.Date("2015-1-1"))
country <- as_Spatial(country[country$gwcode %in% c(475),])

# ... --- load change raster
impr.rs <- readRDS(file.path("data", "nigeria_change_rs.rds"))

# ... plot

# ... --- palette
colpal <- inferno(n = 256, alpha = 1, begin = 0, end = 1, direction = -1)
zlim = c(-22, 12)

# ... --- lagos and abuja
abuja.pt <- SpatialPoints(country@data[,c( "caplong",  "caplat" )])
lagos.pt <- SpatialPoints(matrix(c(3.395833, 6.453056), nrow = 1))


# ... --- plot
stub <- c("reg","nat")
for(i in c(1:2)){
  png(paste0(fig.path, "figa8_",stub[i],"_invar.png"), width = 5, height = 4.25, units = "in", res = 300)
  par(mar = c(0,0,0,0))
  print(raster::plot(impr.rs[[i]], useRaster = T, interpolate = F, axes=FALSE, bty="n", box=FALSE,
                     col = colpal, zlim = zlim, legend = c(T,F)[i]))
  print(points(abuja.pt, col = "black", pch = 19))
  print(points(abuja.pt, col = "white"))
  print(points(abuja.pt, col = "white", pch = 20))
  print(text(abuja.pt@coords, labels = "Abuja", col = "white", pos = 4))
  print(points(lagos.pt, col = "black", pch = 19))
  print(points(lagos.pt, col = "white"))
  print(points(lagos.pt, col = "white", pch = 20))
  print(text(lagos.pt@coords, labels = "Lagos", col = "white", pos = 3))
  dev.off()
}


# FIGURE A9: Plot stacked distributions

# ... load data
incr.dens.df <- read.csv(file.path("data", "statereach_distribution.csv"))


# ... ranges (hard-coded)
ranges.ls <- list(regcap.dist.2016 = c(-5,35),
                  natcap.dist.2016 = c(-5,48),
                  regcap.change = c(-24,24),
                  natcap.change = c(-48,15))

# ... Titles
titles.ls <- list(regcap.dist.2016 = "Time to reg. capital 2016 (hrs) \n  ",
                  natcap.dist.2016 = "Time to nat. capital 2016 (hrs) \n  ",
                  regcap.change = "1966-2016 difference in \n time to reg. capital (hrs)",
                  natcap.change = "1966-2016 difference in \n time to nat. capital (hrs) ")
# ... plot paper
for(t in c("regcap", "natcap")){
  pos.df <- incr.dens.df[incr.dens.df$quantile == 50 &
                           incr.dens.df$cowcode != 1 & 
                           incr.dens.df$type == paste0(t,".dist.2016"), 
                         c("cowcode", "value")]
  pos.df <- pos.df[order(pos.df$value),]
  pos.df$pos <- c(1:nrow(pos.df))
  
  this.incr.dens.df <- join(incr.dens.df[grepl(t, incr.dens.df$type),], 
                            pos.df[,!colnames(pos.df) %in% c("value")], type = "left",
                            by = "cowcode")
  this.incr.dens.df$pos[this.incr.dens.df$cowcode == 1] <- max(this.incr.dens.df$pos, na.rm = T) + 1
  this.incr.dens.df$cowcode <- as.factor(this.incr.dens.df$cowcode)
  # this.incr.dens.df$title <- as.vector(titles.ls)[this.incr.dens.df$type]
  
  
  for(st in c(".dist.2016", ".change")){
    # ... data to plot
    sub.type.df <- this.incr.dens.df[this.incr.dens.df$type == paste0(t, st),]
    
    # ... legend filler
    leg.fill <- data.frame(pos = 0, cowcode = "999", 
                           value = seq(min(sub.type.df$value), 
                                       max(sub.type.df$value), length.out = 100) )
    
    png(paste0(fig.path, "figa9_density_",t, st, ".png"), 
        width = ifelse(st == ".dist.2016",4.5,3.5), 
        height = 10, units = "in", res = 300)
    g <- ggplot(smartbind(sub.type.df,
                          leg.fill), 
                aes(x = value, y = pos, group = cowcode,  fill = ..x..)) + 
      geom_vline(xintercept = 0, col = "grey", lty = 2) +
      geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01,  gradient_lwd =.8, color = "white") +
      scale_fill_viridis(name = "",option = "inferno", direction = ifelse(st == ".change", -1, -1), alpha = 1) + 
      guides(fill=FALSE) + 
      xlim(ranges.ls[[paste0(t, st)]]) + 
      xlab(titles.ls[[paste0(t, st)]]) +
      theme_minimal() +
      NULL
    if(st == ".dist.2016"){
      g <- g + scale_y_continuous(name = NULL,  minor_breaks = NULL,
                                  breaks = unique(sub.type.df$pos), 
                                  labels = unique(sub.type.df$cntry_name), 
                                  limits = c(0, max(sub.type.df$pos) + 2))
    } else {
      g <- g + scale_y_continuous(name = NULL, minor_breaks = NULL,
                         breaks = c(1:max(sub.type.df$pos)), 
                         labels = rep("", max(sub.type.df$pos)), 
                         limits = c(0, max(sub.type.df$pos) + 2)) +
        theme(axis.ticks = element_blank())
      
    }
    print(g)
    dev.off()
  }
}

# ... numbers

# ... --- national capital, median
natcap.median.cowcode <- c(490, 432, 475, 517)
for(c in natcap.median.cowcode){
  writeLines(as.character(round_any(incr.dens.df$value[incr.dens.df$type == "natcap.dist.2016" & 
                                                         incr.dens.df$cowcode == c & 
                                                         incr.dens.df$quantile == 50],.1)) ,
             paste0(num.path,"natcap_median_",c,"_2015.tex"))
}

# ... --- national capital, median, change
natcap.median.cowcode <- c(490, 432, 475, 517)
for(c in natcap.median.cowcode){
  writeLines(as.character(round_any(incr.dens.df$value[incr.dens.df$type == "natcap.dist.2016" & 
                                                         incr.dens.df$cowcode == c & 
                                                         incr.dens.df$quantile == 50],.1)) ,
             paste0(num.path,"natcap_median_diff_",c,"_2015.tex"))
}



# FIGURE A6: Increase in number of admin units ###

# ... load data
adm.cow.yrs <- read.csv(file.path("data", "adminnum_bycountry.csv"))

# ... average by year
adm.agg <- aggregate.data.frame(list(number = adm.cow.yrs$number),
                                list(year = adm.cow.yrs[,c( "year")]), FUN = mean)


# ... plot
png(paste0(fig.path, "figa6_regnum.png"), width = 7/3, height = 7/3, units = "in", res = 300)
g <- ggplot(adm.agg, aes(x = year, y = number)) +
  geom_line(data = adm.cow.yrs, aes(group = cowcode), col = "grey", alpha = .5)  + geom_line() + 
  xlab("Year") + ylab("Number of regions") +
  theme_minimal()
print(g)
dev.off()




# FIGURE A6:  Road Building #######################

# ... load data
built.cow.df <- read.csv(file.path("data", "roadidx_bycountry.csv"))
road.agg <- read.csv(file.path("data", "roadidx_continent.csv"))

# ... plot
png(paste0(fig.path, "figa6_roadsbuilt.png"), width = 7/3, height = 7/3, units = "in", res = 300)
g <- ggplot(road.agg, aes(x = year, y = built)) +
  geom_line(data = built.cow.df, aes(group = COWCODE), col = "grey", alpha = .5) + 
  geom_line() + 
  xlab("Year") + ylab(paste0("Roads, qual.-weighted,\nindexed in 1966")) +
  ylim(.75,2.5) + scale_x_continuous(breaks = pretty(built.cow.df$year, n = 3)) +
  theme_minimal() + 
  NULL
print(g)
dev.off()

# FIGURE A6:  Urbanization ########################

# ... load data
urban.rates <- read.csv(file.path("data", "urbanization_bycountry.csv"))

# ... aggregate to continent
urban.agg <- aggregate.data.frame(urban.rates[,c("SP.POP.TOTL","SP.URB.TOTL")],
                                  list(year = urban.rates[,c("year")]), FUN = sum, na.rm = T)

# ... plot
png(paste0(fig.path, "figa6_urbanpop.png"), width = 7/3, height = 7/3, units = "in", res = 300)
g <- ggplot(urban.agg, aes(x = year, y = SP.URB.TOTL/SP.POP.TOTL)) +
  geom_line(data = urban.rates, aes(group = country.name), col = "grey", alpha = .5)  + geom_line() + 
  xlab("Year") + ylab("Urban population (%)") + scale_x_continuous(breaks = c(1960, 1980, 2000)) +
  theme_minimal() + 
  NULL
print(g)
dev.off()


